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WEIGHTED  LEAST  SQUARES  FIT  OF  A  REAL  TONE  TO  DISCRETE  DATA, 

BY  MEANS  OF  AN  EFFICIENT  FAST  FOURIER  TRANSFORM  SEARCH 

INTRODUCTION 

Estimation  of  the  parameters  of  a  tone  with  unknown  amplitude,  frequency, 
and/or  phase  has  attracted  considerable  attention;  see,  for  example,  [1-9]. 
However,  fitting  data  with  a  single  pure  complex  tone  leads  to  a  simpler 
search  problem  than  fitting  with  a  real  tone  (as  will  be  demonstrated  in  the 
next  section).  In  particular,  fitting  with  a  complex  tone  was  considered  in 
[1-5,  7],  while  fitting  with  real  tones  has  been  the  subject  of  [6,  8,  9], 
However,  the  frequency  of  the  tone  was  assumed  known  in  [6,  8],  whereas  it  had 
to  be  estimated  in  [9]. 

Here  we  will  extend  the  results  in  several  directions  for  the  case  of 
fitting  real  data  with  a  real  tone.  First,  arbitrary  real  weighting  of  the 
errors  at  each  discrete  instant  are  incorporated.  Second,  the  function  that 
must  be  searched  for  a  maximum  is  manipulated  into  a  form  which  requires  that 
only  two  FFTs  of  two  real  sequences  be  conducted.  Third,  these  two  operations 
are  combined  into  one  FFT  of  a  complex  sequence,  the  outputs  of  which  are 
decoupled  in  a  very  efficient  manner,  in  order  to  yield  the  desired  search 
function.  Fourth,  parabolic  interpolation  of  the  three  outputs  in  the 
neighborhood  of  the  search  maximum  is  employed  in  order  to  give  a  refined 
estimate  of  the  tone  frequency.  Finally,  a  minute  search  for  the  best  tone 
frequency  is  conducted,  the  extent  of  which  is  left  up  to  the  user.  The  end 
result  of  this  investigation  is  a  program  for  conducting  an  efficient  and  fast 
fine-grained  search  for  the  determination  of  the  unknown  amplitude,  frequency, 
and  phase  of  the  best-fitting  real  tone  to  a  given  set  of  discrete  real  data 
and  subject  to  any  error  weighting  of  interest. 

This  procedure  is  applicable  to  arbitrary  data  record  lengths.  Also,  no 
assumptions  about  the  statistics  of  any  additive  noise,  that  may  be  present  in 
the  data  record,  are  made.  However,  when  the  available  data  record  is  the 
result  of  a  pure  tone  and  additive  zero-mean  Gaussian  noise,  the  procedure  can 
be  interpreted  as  maximum  likelihood  estimation  [9], 
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ERROR  MINIMIZATION 


Before  we  begin  the  detailed  investigation  of  fitting  a  real  tone  to  real 
data,  we  first  consider  the  simpler  problem  of  fitting  a  pure  complex  tone. 
This  will  serve  as  a  comparison  procedure  and  will  oack  up  the  statement  made 
in  the  Introduction. 

COMPLEX  TONE 

The  discrete  data  available  consist  of  N  values  {xk},  taken  at 
increment  a.  If  the  data  are  complex  and  we  fit  the  data  with  a  pure  complex 
tone,  we  must  address  the  problem  of  minimizing  the  weighted  squared  error 


(1) 


k 


where  the  summation  on  k  is  taken  over  all  nonzero  summands.  Normally,  the 


taken  to  be  nonzero  over  the  range 


1  <_  k  <_  N;  however,  the  presentation  allows  for  any  range  of  the  variable  k. 
The  parameter  a  in  (1)  is  the  complex  amplitude,  and  u>  is  the  pure  tone 
(radian)  frequency,  which  is  presumed  real. 

If  we  consider  u>  given  for  the  moment  in  (1),  the  best  choice  of  a  to 
minimize  error  E  is  given  by 


(2) 


Substitution  of  this  result  for  a  in  (1)  results  in  error 


E(“>  -  2"klxkl2  -  15. 


k 


k 


p(-ioukA) 


k 


(3) 
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This  error  is  minimized  by  choosing  frequency  u>  to  maximize  the  quantity 


(4) 


which  is  the  standard  magnitude-squared  Fourier  transform  of  the  weighted 
data.  Thus,  direct  application  of  an  FFT  is  a  good  procedure  to  apply  to  this 
problem  and  has  been  so  employed  in  the  past  [5].  Since  (4)  has  period  2k/a 
in  uj,  there  is  no  need  to  compute  (4)  except  for  the  range  -u  <_  wA  <_  ir. 

REAL  TONE 

We  now  restrict  consideration  to  the  case  of  major  interest  here,  namely, 
real  data  [x^,  and  attempt  to  fit  it  with  samples  of  a  pure  real  tone,  that 
is, 


a  cos(u)ka)  +  e  sin(uikA)  . 


(5) 


Here,  a  and  6  are  the  real  coefficients  of  the  in-phase  and  quadrature 
components  of  the  tone.  If  we  let  "normalized  frequency" 


a  =  0)A  y 


(6) 


the  weighted  squared  error  to  be  minimized  is 


k 


For  later  use,  we  define  the  two  Fourier  series: 


k 


(8) 


k 
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The  first  is  the  window  associated  with  weights  jw^,  while  the  latter  is 
the  Fourier  transform  of  the  weighted  data. 


The  variable  a  appearing  in  (7)  will  be  called  the  "frequency"  of  the 
tone.  If  we  consider  frequency  a  given  for  the  moment,  setting  the  partial 
derivatives  of  error  E  with  respect  to  a  and  a,  both  equal  to  zero,  results  in 
the  pair  of  simultaneous  linear  equations  for  their  optimum  values: 


A11  °o  +  A12  eo  =  Lr(a)  • 


A12  “o  +  A22  8o  “  Li ^ 


(9) 


Here  sub  r  and  i  denote  real  and  imaginary  parts,  respectively.  We  also  have 
the  scale  factors  expressible  in  the  forms 


A11  =  2wk  cos2(ak)  =  7&K0)  +  wr(2a)~]  » 

k 

A22  =  sin2(ak)  =  ^[w(0)  -  Hr(2a)J  , 

k 

A^2  =  cos(ak)  sin(ak)  =  *“  ‘g’  Wi  C2a)  ,  (10) 

k 

where  we  have  made  extensive  use  of  (8).  Solution  of  (9)  yields  for  the  tone 
coefficients. 


ao 


A22Lr^  +  A12Li  ^ 

- _ - 

A  A  -  A 
H11H22  rt12 


6o  " 


-AnL.(a)  -  A12Lr(a) 
A11A22  "  A12 


(ID 


The  use  of  ( 9 )— ( 11 )  in  error  (7)  now  results  in  modified  error 
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E ( a )  *  ^wk  [xk  -  aQ  cos(ak)  -  bq  sin(ak)]2  = 
k 

=  2wk  ^-xk  "  ao  cos(ak)  “  s0  sin(ak) ]  xk  = 
k 


-  xk  -  “o  Lr(a>  +  6o  Li<a>  “ 

k 

■  Zwk  xk  -  B<a>  • 

k 


where  we  define  real  function 


(12) 


B(a)  =  aQ  Lr(a)  -  eQ  L^a)  = 

A22L^(a)  +  2A12Lr(a)L.  (a)  +  AuL2(a) 

A  A  _ 

11*22  *12 

This  quantity,  which  must  now  be  maximized  by  choice  of  a,  was  previously 
encountered  in  [9;  (10)],  but  limited  there  to  the  case  of  equal  weights 
{wkl  We  concentrate  henceforth  on  function  B(a),  aware  that  we  can  always 
return  to  error  E(a)  by  means  of  (12). 
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MANIPULATIONS  OF  8(a) 


In  this  section,  we  derive  alternative  forms,  properties,  and 
interpretations  of  the  function  B(a).  Tne  weighted  squared  error  is  directly 
related  to  B(a)  by  means  of  (12). 

ALTERNATIVE  FORM  FOR  B(a) 


A  more  useful  and  compact  form  for  B(a) 
to  (10)  reveals  that  the  denominator  of  (13) 


in  (13)  is  possible, 
is  simply 


Reference 


^[w2(0)  -  |w2(2a)|]  .  (14) 

Similarly,  use  of  (10)  allows  development  of  the  numerator  of  (13)  according  to 


|[w(0)  -  Wr(2a)]L2(a)  -  W.  (2a)Lr(a)Li  (a)  +  |[w(0)  +  Wr(2aj]  L2(a)  = 

=  \  W(0)|L2(a)|  -  ^[wr(2a)L2(a)  ♦  2Wi  (2a)Lr(a)Li  (a)  -  Wr(2a)L^(a)]  = 

=  \  W(0)  |l 2 ( a ) |  -  \  Re{w*(2a)L2(a)}  .  (15) 

Coupling  (14)  and  (15)  together,  the  expression  in  (13)  becomes 

6(a)  --  2  W(0)lL2U)l-Re{/(2a)L2(a)}  _  (16) 

w  (0)  -  |w^(2a)| 

The  required  quantities  here  are  available  from  (8)  as 

W(2a)  =  K  exp(-i2ak), 
k 

L(a)  =  ^>wk  xk  exp(-iak)  .  (17) 

k 
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It  is  immediately  obvious  from  (16)  that  8(a)  can  never  be  negative 
(presuming  that  the  weights  are  nonnegative). 

The  general  result  for  8(a)  in  (16)  is  the  quantity  that  must  be 
maximized  by  choice  of  frequency  a.  However,  it  is  interesting  to  observe 
that  for  frequencies  where  the  window  is  small,  that  is, 

|w(2a)|  «  W (0 )  ,  (18) 


then  (16)  simplifies  to 


B(a)  z 


xk  exp(-iak) 


(19) 


which  is  identical  to  (4).  Thus,  for  those  frequencies  where  (18)  is  true, 
the  function  B(a)  is  approximately  the  magnitude-squared  Fourier  transform  of 
the  weighted  data;  this  corresponds  to  values  of  a  not  near  multiples  of  w. 


PROPERTIES  OF  B(a) 


Since  W(2a)  has  period  n  in  a,  while  L(a)  has  period  2ir  in  a,  the 

/ 

function  B(a)  in  (16)  must  have  period  2tt  in  a;  that  is, 

B(a  +  2tt)  =  8(a)  .  (20) 

But  at  the  same  time,  we  have  even  property 

B(-a)  =  B(a)  ,  (21) 

★  ★ 

because  L(-a)  =  L  (a),  W(-2a)  =  W  (2a),  using  the  realness  of  sequences 
{wk}  and  fxk}.  What  this  means  is  that  we  only  need  to 

compute  B(a)  for  0  <  a  £  u  ,  (22) 

since  all  other  values  can  be  obtained  therefrom.  Reference  to  (6)  reveals 
that  u  is  being  varied  over  the  range  (0,  -n- / a ) ,  or  that  cyclic  frequency 
f  =  to/ ( 2 tt )  is  varying  over  (0,  .5/a).  This  latter  range  extends  up  to  the 
Nyquist  frequency,  as  expected. 
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VALUE  OF  B ( a )  AS  a  »  0 

If  we  substitute  a  =  0  in  (16),  we  get  B(0)  =  0/0,  which  is 
indeterminate.  Hence,  for  small  a,  we  expand  (17)  according  to 

W(2a)  ~  Wq  -  i2aW^  -  2a2W^  , 


1  2 

L(a)  ~  Lq  -  iaLj^  -  a  L2  , 


(23) 


where  n-th  order  "moments" 


W  =  Sw.  kn  , 

n  <C  k  ’ 

k 

Ln  '  £wk  xk  k"  ’ 
k 

Substitution  of  (23)  in  (16)  and  simplification  yields 


(24) 


lim  B(a) 
a»0 


vi  -  2w0; 

U0W2  ”  “l 


This  limiting  result  is  the  same  value  that  is  attained  as  if  we 
minimized  weighted  error 


(25) 


E  =  ^  wk  [xk  -  U  -  vk]2  ,  (26) 

k 

by  choice  of  constant  value  u  and  linear  trend  vk.  In  fact,  direct 
minimization  of  (26)  yields  optimum  coefficients 


-  W2L0  W1L1 

^  ^ 

W0W2  1 


v0  = 


W0L1  ~  Wl4) 

V2  -  “l 


and  associated  minimum  error 


(27) 
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=  1" 


k  xk 


W  L2 

0  1 


2W1L1L0  +  W2L0 


WqW2  -  Wj 


(28) 


As  claimed  above,  the  last  term  in  (28)  is  precisely  the  result  given  by  (25); 
see  (12)  also.  Thus,  the  limit,  as  a  >  0,  of  model  fit  (7)  is  the  best¬ 
fitting  constant  plus  linear  trend  to  the  given  data.  This  can  be  obtained 
from  (7)  only  if  quadrature  coefficient  e  behaves  as  1/a  as  a  >  0.  Indeed,  in 
a  later  section,  we  will  show  that  this  is  precisely  the  behavior  of  6  in  this 
limit.  Thus,  setting  a  =  0  in  (7)  and  keeping  &  finite  does  not  lead  to  the 
result  in  (25)  and  (28),  but  instead  gives  only  the  best  fitting  constant.  We 
will  allow  the  more  general  fit  afforded  by  (26)  here,  and  will  utilize  the 
value  achieved  by  (25)  in  the  limit,  as  a  »  0. 


VALUE  OF  B(a)  AS  a  »  it 


If  we  substitute  a  =  it  in  (16),  there  follows  B(ir)  =  0/0,  which  is 
indeterminate.  However,  if  we  let  a  =  ir  +  a,  we  see  from  (17)  that 


W(2a)  =  W(2tt  +  2a)  =  W(2a)  , 

L(a)  =  L(ir  +  a)  =  wk  (-l)k  xk  exp(-iak)  .  (29) 

k 

Thus,  W(2a)  behaves  the  same  about  a  =  0  as  W(2a)  does  about  a  =  0.  Also,  the 
last  term  in  (29)  behaves  the  same  about  S’  =  0  as  L(a)  does  about  a  =  0, 

Jy 

provided  that  each  data  element  xk  is  replaced  by  (-1)  xk.  Thus,  (25) 
can  be  immediately  utilized  to  yield  the  result 


1 im  B(a) 

a>7r 


“oLl  -  2“l4[0  * 


“0U2  -  “l 


where  moments  (24)  have  been  replaced  by 


(30) 


(31) 
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Physical  interpretation  of  result  (30)  is  similar  to  that  given  earlier 
for  a  >  0  in  ( 25 ) — ( 28) .  Namely,  in  the  limit  as  a  >  it,  the  best  constant  plus 
linear  trend  is  fitted  to  alternating  data  ^(-lj  x^"].  Again,  this  requires 

quadrature  coefficient  e  in  fit  (7)  to  behave  like  1  / ( a  -  ir)  as  a  >  n. 

EXAMPLE  OF  EQUAL  WEIGHTS 

Let  weights 


(32) 


Then  window  (17)  becomes 


exp(-i(N  +  l)a)  . 


(33) 


This  is  the  example  considered  in  [9]. 


The  moments  (24)  for  this  case  are  given  by 


Wq  -  1,  -  |(N  +  1),  W 2  =  £(N  +  1)(2N  +  1)  , 


WQW2  -  W2  .  ^(N2  -  1)  . 


(34) 


The  numerator  of  (25)  is  then 


L^  -  (N  +  IJLjLq  +  |(N  +  1) (2N  +  1)L2  = 


(35) 
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Then  (12),  (25),  and  (34)  yield 


1  im 
a>0 


1  52  i_ 

r^xT  12 

r-c  fk  n  +  iX] 

N  f  k  N2 

k 

kJ  N2(N2  -  1) 

Lf  k(  '  2  7J 

(36) 


which  can  be  recognized  as  the  minimum  error  for  the  best-fitting  constant 
plus  linear  trend  to  data  . 
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IN-PHASE  AND  QUADRATURE  COEFFICIENTS 

Since  the  modeling  waveform  in  (7)  is 

a  cos(ak)  +  s  sin(ak)  =  Re  {(a  -  is)  exp(iak)}  , 

the  complex  coefficient  or  strength  of  pore  complex  tone  exp(iak)  is  a 
From  (11)  and  (10),  the  numerator  of  aQ  -  isQ  is  expressible  as 

A22Lr^  +  A12Li^  +  iAllLi(a)  +  iAi2Ma)  = 

=  ?[w(0)  -  Wr(2a)]Lr(a)  +  i|[w(0)  +  Wr(2aj]L.(a)  + 

+  1  ("  |]Wi(2a)[Lr(a)  -  iL. (a)~j  = 

=  ~  W(0)L(a)  -  \  W(2a)L*(a)  .  (38) 

Combining  this  with  the  denominator  previously  computed  in  (14),  we  have  for 
the  optimum  complex  coefficient, 

?  W(0)L(a)  -  W(2a)L*(a)  (39) 

0  0  W^(0)  -  |vr(2a)| 

For  frequencies  a  such  that  the  window  is  small  relative  to  the  origin 
value  (see  (18)),  (39)  simplifies  to  the  approximate  result 

aQ  -  iB0S  2L(a)  =  2  £wk  xk  exp(-iak)  ,  (40) 

k 

which  is  just  the  Fourier  transform  of  the  weighted  data. 

If  only  the  phase  of  the  real  tone  is  of  interest,  (39)  indicates  that 


(37) 
-  is. 
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arg(a0  -  ieQ)  =  arg(w(0)L(a)  -  W(2a)L*(a))  .  (41) 

If  frequency  a  is  known,  this  result  is  directly  applicable;  but  if  a  is 
unknown,  the  value  a  that  maximizes  (16)  must  be  used. 

NORMALIZATION  OF  WEIGHTS 

Without  loss  of  generality,  the  sum  of  the  weights  ^w^  can  be  set 
equal  to  unity;  that  is,  set 


W(0)  =  5  wk  =  1  . 

k 

Then  the  complex  coefficient  in  (39)  reduces  to 


a 


0 


.  =  2  L( a)  -  W(2a)L*(a) 

°  1  -  |w2(2a)| 


while  the  maximizing  function  B(a)  in  (16)  becomes 


(42) 


(43) 


B(a)  =  2  lL^(a)l  -  Re  { W*(2a)L2(a j) 
1  -  |w2(2a)| 


(44) 


This  slightly  reduces  the  number  of  computations  that  have  to  be  conducted  and 
has  been  adopted  in  the  program  written  here.  This  scaling  is  also  retained 
in  the  following  subsection. 


INTERPRETATION  OF  (43) 

An  alternative  form  for  coefficient  (43)  is 

.  o  L(a)  -  W(2a)L(-a) 

»  i  -  plilli  1 


(45) 
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where  we  utilized  the  realness  of  data  and  weights  £wkj  .  This  result 
can  be  interpreted  as  follows:  the  term 


2L(a)  =  2  ^wk  \  exp(-iak) 
k 


(46) 


is  an  estimate  of  the  complex  strength  of  the  positive-frequency  complex 
exponential  exp(iak)  in  the  real  data  {xk}  ,  as  modified  by  the  weights. 
Similarly,  2L(-a)  estimates  the  strength  of  the  term  exp(-iak)  in  the  real 
data.  The  window  W(2a)  measures  the  amount  of  spillover  from  frequency  -a  to 
frequency  a,  that  is,  at  separation  2a,  due  to  the  weights  [w^.  This 
fraction  (including  phase  information)  of  the  spillover  from  negative 
frequencies  to  positive  frequencies  is  subtracted  from  strength  2L(a). 
Finally,  the  denominator  factor  1  -  |w2(2a)|  renormalizes  the  remainder 
according  to  the  fractional  spillover. 

To  justify  this  last  scale  factor,  suppose  that  the  data  (xk"J  contain  a 
pure  real  tone  at  precisely  the  frequency  a;  that  is,  let 


x,  =  a  cos(ak)  +  sin(ak)  = 
k  o  o  '  ' 


=  Re {(aQ  -  ieQ)  exp(iak)}  . 


(47) 


Then  (46)  yields 


2L(a)  =  ^  wk[a0(exp(iak)  +  exp(-iak))  -  ieo(exp(iak)  -  exp(-iak)  )]exp(-iak)  = 
k 


a0[l  +  W(2a)]  -  ieQ[l  -  W(2a) ]  . 


(48) 


Therefore, 


2L(-a)  =  2L*(a)  =  a0[l  +  W*(2a)]  +  i SQ[1  -  W*(2a)]  . 


(49) 


Therefore,  the  numerator  of  aQ  -  i6Q  in  (45)  is 
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2L (a)  -  W(2a)2L(-a)  = 

=  aQ(l  +  W)  -  i 0O(1  -  W)  -  W[a0(l  +  W*)  +  iS0(l  -  W*)]  = 

=  ^ao  “  i0o)(1  “  lw2|)  =  K  "  Uo^1  “  lw2(2a)|  3  ,  (50) 

where  we  adopted  the  notational  simplification  W  =  W(2a)  during  the 
manipulations.  Thus,  the  denominator  factor  1  -  Jw^(2a)|  in  (45)  is 
necessary  to  scale  the  amplitude  back  up  to  its  correct  value  of  aQ  -  ieQ. 

VALUE  OF  COEFFICIENT  AS  a  >  0 


We  want  to  investigate  the  behavior  of  coefficient  aQ  -  ieQ  in  (39) 
as  a  »  0.  (If  we  try  to  set  a  =  0,  we  get  <*0  -  i60  ■  0/0,  which  is 
indeterminate.)  Accordingly,  substitute  expansions  (23)— ( 24)  into  (39)  and 
simplify  to  obtain  the  expression 


ao  '  10d^ 


W2l0  Wl4 

«0W2  - 


i  woLi  -  W1L0 

S  Wz  -  “! 


as  a  »  0 


(51) 


This  result  corroborates  the  claim  made  under  (28)  that  eQ  behaves  as  1/a  as 
a  >  0.  That  is,  the  optimum  quadrature  coefficient  of  the  pure  real  tone  gets 
arbitrarily  large  as  frequency  a  tends  to  zero. 


If  we  combine  (51)  with  the  modelling  function  in  (7),  we  have 


ctQ  cos(ak)  +  &o  sin(ak)  ~  aQ  +  &Qak  ~ 


W  L  -  W  L 
20  W1L1  , 

- 2“ 

W0W2  “  W1 


W0L1  Ul4) 

W0W2  "  W1 


as  a  »  0 


(52) 


which  is  precisely  (26)  and  (27).  Thus,  the  limit,  as  a  >  0,  of  modeling  (7) 
is  to  fit  the  best  constant  plus  linear  trend  to  the  data. 


15 


TR  7785 


FF r  REALIZATION 


For  purposes  of  minimizing  computations,  we  henceforth  assume  that  the 
weights  have  been  normalized  according  to  (42);  that  is,  their  sum  equals 
unity.  This  feature  is  incorporated  in  the  following  equations  and  the 
resultant  program. 

MANIPULATION  INTO  FFT  FORMS 

According  to  (22),  we  are  interested  in  evaluating  B(a)  in  (44)  over  the 
range  0  £  a  £  it,  where  functions  W  and  L  are  given  by  (17).  Suppose  then  that 
we  focus  attention  on  values  of  frequency  a  given  by 


(53) 


Integer  M  will  be  chosen  to  be  a  power  of  2,  and  is  unrelated  to  N,  the  number 
of  data  points.  Then  (17)  yields 


(54) 


k 


which  is  recognized  as  an  M-size  FFT  of  N  nonzero  real  weighted  data  values 


At  the  same  time,  (17)  also  gives  for  the  window 


1 


vi. 12  exp(-i2nmj/M)  , 


(55) 


j  even 


where  we  let  j  =  2k.  Now  if  we  define  sequence 


16 


d.  = 
J 


wj/2  for  j  even 


for  j  odd 


then  (55)  becomes 


W  =  exP(-i2irmj/M)E  J  (m)  , 

j 

which  is  an  M-size  FFT  of  sequence  {d.|. 

J 

Direct  employment  of  (54)  and  (57)  in  (44)  yields 
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(56) 


(57) 


U2(m) 

|  -  Re  ' 

>*(■")  z2^)! 

1  - 

1 

(58) 


Thus,  if  we  evaluate  the  two  FFTs  for  (<£(m)]  and  £$(m)}  in  (54)  and  (57), 
respectively,  we  have  all  the  quantities  necessary  to  determine  B(m2ir/M)  for 
0  <  m  <  M/2. 


TWO  REAL  FFT's  VIA  ONE  COMPLEX  FFT 


Since  (54)  and  (57)  constitute  FFTs  of  real  sequences,  they  are  not 
making  full  use  of  the  capabilities  of  an  FFT.  To  exploit  the  inherently 
complex  nature  of  this  tool,  let 

Zk  =  wkxk  +  idk  for  1  <  k  <  2N  ,  (59) 

where  sequence  £d^  was  defined  in  (56).  (Half  of  the  real  terms  and  half 

of  the  imaginary  terms  are  zero  in  (59).)  Then  the  FFT  of  size  M  of  sequence 
(59)  is 

Z(m)  =  ^zk  exp(-i2irmk/M)  =  <£(m)  +  i  3Tm)  ,  (60) 

k 

where  we  presume  that  M  >  2N.  (Methods  of  circumventing  this  limitation  are 
given  in  [10].) 
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Using  the  realness  of  sequences  [x^j,  ^  f°^ows  from 

(54)  and  (57)  that 


Z  (M  -  m)  =  (M  -  m)  -  i#  (M  -  m)  =  ^(m)  -  ij)(m)  .  (61) 

Now  combining  (60)  and  (61),  we  have 

2  £(m)  =  Z(m)  +  Z  (M  -  m)  =  Sx(m)  +  i  Dy(m)  , 

2a&(m)  =  -i[Z(m)  -  Z  (M  -  in)]  =  Sy(m)  -  i  Dx(m)  ,  (62) 

where  the  real  sum  and  difference  functions  are  defined  as 

Sx(m)  =  X(m)  +  X(M  -  m), 

S  (m)  =  Y(m)  +  Y(M  -  m), 

Dx(m)  =  X(m)  -  X(M  -  m), 

Dy(m)  =  Y (m)  -  Y(M  -  m),  (63) 

in  terms  of  the  real  and  imaginary  parts  of  FFT  output  Z(m)  in  (60),  namely. 


Z(m)  =  X(m)  +  i  Y(m)  . 


(64) 


Equation  (62)  accomplishes  the  decoupling  of  the  FFT  output  Z(m)  so  as  to 
yield  the  two  desired  FFTs  i£(m)  and  £(m)  indicated  in  (54)  and  (57). 

However,  it  is  advantageous  to  continue  with  the  breakdown  of  these  two 
complex  sequences  <i£(m)  and  ^(m),  as  done  in  (62),  in  terms  of  all  the  purely 
real  quantities  given  in  (63).  For  upon  substitution  of  (62)  in  desired 
quantity  B(m2u/M)  in  (58),  we  obtain  the  simplified  form 


B 


Sx(m)[2  -  S  (m)]  +  D^(m)[2  +  S  (m)]  +  2Sx(m)Dx(m)D  (m) 
- 4  -  $(.)  >  #.)) - 


.  (65) 


This  latter  form,  which  utilizes  only  real  arithmetic,  can  be  used  only  for 
0  <  m  <  M/2.  The  values  for  B(0+)  and  B(n-)  must  come  from  (25)  and  (30), 
respectively,  with  Wg  =  1.  A  program  for  calculation  of  fB(m2ir/M)}  by  means 
of  (56),  (59),  (60),  (63),  and  (65)  is  furnished  in  the  appendix. 
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SELECTION  OF  FFT  SIZE  M 

It  was  presumed  in  (59)  and  (60)  tnat  FFT  size  M  >  2N,  where  N  is  the 
number  of  data  points,  in  order  that  there  be  an  array  element  in  location  M-l 
available  to  receive  data  element  id^.  However,  there  is  an  additional 
reason  for  choosing  M  this  large,  having  to  do  with  the  rate  at  which  B(a) 
varies.  The  function  B(a)  in  (44)  depends  critically  on  window  function 
N(2a).  For  equal  weights,  the  results  in  (32)  and  (33)  indicate  that  W(2a) 
changes  significantly  in  an  interval  of  length  ir/N;  in  fact,  this  is  the 
separation  between  zero  crossings.  If  order  to  track  this  rapid  variation  in 
W,  the  increment  2u/M  in  frequency  a  in  (53)  and  (58)  must  be  smaller  than 
tt / N .  Thus,  requirement  M  >  2N  is  a  minimal  requirement;  in  fact,  it  may  be 
advantageous  to  consider  M  several  times  larger  than  N,  if  storage  and  FFT 
execution  time  are  not  excessive.  Of  course,  the  larger  M  is  taken,  the  less 
fine-grain  interpolation  will  be  required  later. 

For  other  weightings  than  flat,  such  as  Hanning,  where  the  effective 
length  of  the  weighting  is  foreshortened  due  to  taper  at  the  edges,  the  window 
function  W  is  broader,  and  the  condition  on  M  is  alleviated  somewhat. 

However,  M  >  2N  is  a  good  rule  of  thumb  to  use  in  most  cases. 

INTERPOLATION  PROCEDURE 

When  the  complete  set  of  values  of  B(m2u/M)  for  0  £  m  <_  M/2  are 
available,  they  are  searched  to  find  the  maximum  value.  This  maximum  value 
and  the  two  neighboring  bin  outputs  (m  values)  are  then  used  in  a  parabolic 
interpolation  procedure  to  refine  the  estimate  of  the  location  of  the  best 
value  of  frequency  a  and  the  corresponding  maximum  value  of  B(a). 

Finally,  this  latter  value  of  a  can  be  used  as  a  starting  value  for  a 
fine-grained  search,  again  by  means  of  parabolic  interpolation,  in  the 
neighborhood  of  this  peak.  These  features  are  all  incorporated  in  the 
accompanying  program  for  this  search  procedure,  where  direct  use  of  (44)  is 
made;  the  previous  FFT  results  are  of  no  use  in  this  final  vernier 
estimation.  Along  with  each  estimated  frequency  a,  the  corresponding 


19 


TR  7735 


coefficient 
the  vernier 
accuracy  of 


a  -  is  is  also  estimated  and  printed  out.  A  few  stages  of 
oo 

analysis  suffice  to  give  stable  frequency  estimates  within  the 
the  computer  used  here. 
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RESULTS 


An  example  of  N  =  25  data  points  with  FFT  size  M  =  1024  is  displayed  in 
figure  1,  for  the  data  sequence 


xk  =  cos(k)  +  j  sin(k)  for  1  <  k  <_  N 


(66) 


and  for  the  equal  (or  flat)  weighting  case  of  (32).  The  abscissa  is 
normalized  frequency  a  =  wA,  and  the  ordinate  is  B(a)  normalized  relative  to 
its  peak  value.  The  low-level  sidelobes  in  figure  1  are  due  to  the  nature  of 
the  window  W(2a),  given  by  (33)  for  this  case. 

The  line  labeled  INITIAL  gives  the  bin  number  Js  in  which  the  peak  is 
located.  This  bin  and  the  two  adjacent  ones  are  then  interpolated  by  means  of 
a  parabola  to  yield  the  initial  value  for  B(a)  labeled  Big  and  the  abscissa 
estimate  a  =  1.0000445.  This  value  of  a  is  then  employed  in  subroutine  SUB  B 
to  give  the  corresponding  value  B(a). 

In  the  next  two  lines  of  the  print  out,  the  above  value  of  a  is  perturbed 
by  ±  Delta,  the  function  B(a)  is  computed,  and  parabolic  interpolation  is 
again  used  on  these  three  points  to  give  the  estimates  labeled  as  REFINED 
values.  Then  this  refined  a  value  is  used  to  recompute  B(a)  and  indicated  as 
the  MAXIMUM  value  in  the  print  out.  Finally,  the  coefficient  estimates  a  and 
S,  along  with  the  minimum  error,  E  .  ,  are  printed  out. 

The  whole  cycle  of  perturbation  and  parabolic  interpolation  is  repeated 
in  the  next  separated  four  lines  of  print  out,  but  this  time  with  Delta 
decreased  by  a  factor  of  10.  This  cycle  is  repeated  one  final  time  in  each  of 
the  figures  presented.  Prolonged  repetition  would  result  in  excessive 
round-off  error,  due  to  the  differencing  of  similar  function  values. 

If  the  weighting  is  changed  to  Hanning, 


for  1  <_  k  <_  N  , 


(67) 
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the  corresponding  results  are  displayed  in  figure  2.  The  initial  estimate  is 
a  =  .99999970,  which  is  then  refined  to  a  =  1.  The  coefficients  converge 
rapidly  to  the  correct  values,  and  the  figure  displays  no  visible  sidelobes 
for  this  case  of  Hanning  weighting.  However,  the  window  is  broader. 

The  results  of  figures  3  and  4  correspond  to  figures  1  and  2, 
respectively,  except  that  white  noise  of  power  1/12  has  been  added  to  the 
waveform  of  (66).  Now  the  Hanning  weighting  result  in  figure  4  also  displays 
sidelobes,  due  to  random  fitting  of  the  particular  noise  samples  utilized. 

The  refined  values  of  a  converge  to  a  =  .99318  and  a  =  1.00208,  respectively, 
which  are  not  exactly  correct,  due  to  the  additive  noise.  Also,  the 
coefficient  a  and  &  are  considerably  off  their  correct  values,  although  the 
Hanning  results  in  figure  4  are  better  than  for  the  flat  weighting  used  in 
figure  3. 

Figures  5  and  6  are  conducted  for  the  two-tone  data  sequence 

x^  =  cos(k)  +  j  sin(k)  +  cos(2k)  ,  (68) 

with  no  additive  noise.  The  second  peak  near  a  =  2  in  these  figures  is  due  to 
the  attempted  match  of  model  (7)  to  the  data,  when  a  is  near  2.  The  program 
locks  onto  the  stronger  tone  and  indicates  its  frequency  as  a  =  .99447  and 
a  =  1.00128,  respectively.  Again,  the  estimates  of  frequency  a  and 
coefficients  a  and  e  are  better  for  the  Hanning  weighting  in  figure  6  than  for 
the  flat  weighting  in  figure  5.  This  is  due  to  the  lower  sidelobes  of  window 
W(2a)  in  the  Hanning  case. 
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NUMBER  OF  DATA  POINTS  N  =  25 
SIZE  OF  FFT  M  =  1024 

INITIAL:  Js  =  163  Big  =  . 62 1 874542623  a  =  1.00004447034  B<a>  =  . 62107486 


REFINED 

a  — 

1 . 00000000426 

REFINED 

B  <  a  > 

=  . 621074930024 

MAXIMUM 

B(a) 

=  .621874930049 

A  1  p  h  a  = 

.  9  9  9  9  9  9  9  6  9  6  8  5  B  e  t  a  = 

.  5  8  8  0  8  8  8  5  4  3  0  3 

Emin  =  5 . 55 1 1 1 5 1 23 1 3E- 1 6 

REFINED 

a  — 

1 . 0000000001 1 

REFINED 

E  <  a) 

=  . 621074930049 

MAXIMUM 

E  <  a  > 

=  .621074930049 

A  1  p  h  a  = 

.  9  9  9  9  9  9  9  9  9  1  9  9  B  *=  t  a  = 

.  5  0  0  0  0  0  0  0 1436 

Emin  =  1 .  1 1 022302463E- 1 6 

REFINED 

a  - 

1 

REFINED 

B  '•  a  > 

=  .621074930049 

MAXIMUM 

B  ( a  > 

=  .621074930049 

A  1  p  h  a  = 

1 

Beta  =  .49999999  9 

99  4  E  m i n  =  1 . 

,  1 1 0  2  2  3  0  2  4  6  3  E  -  1  6 

Figure  1.  Flat  Weighting,  No  Noise 
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NUMBER  OF  DATA  POINTS  H  =  25 
SIZE  OF  FFT  M  =  1024 
INITIAL:  Js  =  163  Big  =  .62475 
REFINED  a  =  .999999999908 
REFINED  E  < a >  =  .624753873395 
MAXIMUM  B ( a )  =  .624753873395 
Alpha  =  1.08000000059  Beta 

REFINED  a  =  .999999999999 
REFINED  BCa)  =  .624753873395 
MAXIMUM  BCa)  =  .624753873395 
A  1 p h a  =  1.0 8  8  8  0 0 0 0 0 0 1  B  e t a 


REFINED  a  =  .999999999997 
REFINED  BCa)  =  .624753373395 
MAXIMUM  BCa)  =  .624753873395 
Alpha  =  1.00000000802  Beta 


4556  a  =  .999999708554  BCa)  =  .62475387 


.49999999381  Emin  =  8 


.499999999986  Emin  =  -l . 1 1 822382463E- 1 6 
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Figure  2 . 
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HUMBER  OF  BflTfl  POINTS  H 
SIZE  OF  FFT  M  =  1024 


INITIRL: 

J  s  = 

1 62  Big  =  .657123186854  a  =  .993 

221368141 

B < a )  =  .65712576 

REFINED 

a  - 

993 1 78780656 

REFINED 

B  <  a) 

=  .657125332743 

MAXIMUM 

E  a  > 

-=  .  657125832 7 76 

A  1  p  h  a  — 

1.11 

93876303  Beta  =  . 232693423946 

E  rn  i  n  =  6 

. 3994944  8 3 1 4  E - 0  2 

REFINED 

a  = 

99317377595 

REFINED 

B  <  a  > 

=  .657125332776 

MAXIMUM 

E  <  a  > 

=  .657125832776 

A 1  p  h  a  = 

1.11 

9  3 8765  0 0 5  Beta  =  .2  8 2 6  0 8 3  63445 

E  m  i  n  = 

6 . 399494483 1 4E-02 

REFINED 

a  —  . 

993 1 7877584 

REFINED 

B  <  a  > 

=  .657125832776 

MAXIMUM 

B  ( a  ) 

=  . 657 125332776 

A  1  p  h  a  = 

1.11 

9  3  8765 0 51  Beta  =  .2326  0 8 36191 1 

E  rn  i  n  = 

6 . 399494488 1 4E-02 

Figure  3 . 
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Normalized  Frequency  a 


NUMBER  OF  DATA  POINTS  N  =  25 


SIZE  OF 

FFT  M  =  1024 

INITIAL; 

Js  =  1  S3  Big  =  . 64394 

3182613  a  =  1 . 8G 

1 2  0 8  436495  E < 

a.>  =  .64  3  9  4  3  1 

REFINED 

a  =  1 . 80208395996 

REFINED 

B  <  a  >  =  . 643943 178345 

MAXIMUM 

B  < a  >  =  .6 43943178345 

A  1  p  h  a  = 

1.0636516553  Be t a  - 

. 396291769 3 4  2 

E  m in  =  5.15 

5548 1 3299 E -02 

REFINED 

a  =  1.00208395965 

REFINED 

B { a  >  =  .643943178345 

MAXIMUM 

B  < a  >  =  .643943178345 

A  1  pi  h  a  = 

1 . 06365165 68  5  B  e t a 

=  . 396291765184 

Emin  =5.1 

55548 13299E-0 

REFINED 

a  =  1.00208395966 

REFINED 

B  < a  >  =  .643943178345 

MAXIMUM 

B a )  =  .643943178345 

A  1  f  j  h  a  = 

1 . 0  6  3  6  5 1 6  5  6  8  4  B  e  %  a 

=  . 39629176521 1 

Em i n  =  5 . 1 

55548 1 3299E-8 

Figure  4.  Hann 

ing  Weighting 

,  Rdditive 

No  i  se 
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Normalized  Frequency  a 


NUMBER  OF  HfiTfi  POINTS  N  =  25 
SIZE  OF  FFT  M  =  1024 


INITIRLs 

J  s  =  162  Big  =  .6123071 

34756  a  =  . 9945 

19907011 

B <  a)  =  .612 3 056 

REFINED 

REFINED 

MAXIMUM 

A  1  p  h  dL  ~ 

a  =  .994474013322 

B  <  a  >  =  .612305684026 

B(a)  =  .612305684054 
1.03244095052  Beta  = 

. 417726698818 

E  rn  i  n  = 

. 499913742067 

REFINED 
REFINED 
MAX  I  MUM 
Alpha  = 

a  =  .994474009397 

B  C a >  =  . 612305684054 

B < a  >  =  .612305684 0 5 4 

1 . 03244097595  Beta  = 

. 417726642805 

E  rn  i  n  = 

. 499913742 0 6 7 

REFINED 

REFINED 

MAXIMUM 

A  1  p  h  a  — 

a  =  .994474009277 

E  c.  a )  =  .  612305684054 

B<a>  =  .612305684054 
1.03244097673  Beta  = 

. 417726641086 

E  r  n  i  n  = 

. 4999 1 3742067 

Figure  5.  Flat  Weighting,  Two  Tones 
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NUMBER  OF  DATA  POINTS  H  =  25 
SIZE  OF  FFT  M  =  1824 

INITIALS  Js  =  163  Big  =  .623356825941 
REFINED  a  =  1.88127551942 
REFINED  B  <  a  >  =  . 623356862532 
MAXIMUM  B  (  a)  =  .623356862582 
Alpha  =  .989782416869  Beta  = 


1 . 88  1 


51724932123 


REFINED  a  =  1.88127551953 
REFINED  B < a)  =  . 623356862582 
MAXIMUM  B ( a )  =  . 623356862582 
A  1 p h a  -  .9897  8 2415852  B  e  t a  = 


.51724932317 


REFINED  a  =  1.88127551958 
REFINED  B < a )  =  . 623356862582 
MAXIMUM  E  < a  >  =  .623356862532 
Alpha  =  .989782415838  Beta 


.517249323284 


7468591  B(a>  =  .623: 


E  m in  =  .49  9  95856945  8 


Emin  =  .49995856945 


E  m in  =  .49  9  9  5  8  56945 


56 


Figure  6.  Hanning  Weighting,  Two  Tones 
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SUMMARY 

An  automatic  procedure  for  determining  the  best  frequency,  amplitude,  and 
phase  of  a  real  tone  fitted  to  discrete  real  data  has  been  devised  and 
programmed.  It  employs  a  single  complex  FfT  for  the  initial  search  and  then 
refines  the  estimates  by  simple  parabolic  interpolation  procedures.  The  size 
M  of  the  FFT  is  unrelated  to  the  number  N  of  data  points,  but  should  generally 
be  taken  at  least  equal  to  2N  in  order  to  guarantee  adequate  sampling  in  the 
frequency  search.  The  user  can  input  any  real  weights  of  his  choosing 

into  the  program;  these  are  then  automatically  normalized  to  make  their  sum 
equal  to  unity. 

The  procedure  is  applicable  to  data  records  of  any  length  N,  without  any 
approximations.  However,  if  there  is  considerable  noise  in  the  data,  then 
large  N  will  be  required  in  order  to  attain  accurate  estimates  of  the  tone 
frequency,  amplitude,  and  phase.  This  is  not  a  drawback  of  the  least  squares 
procedure  or  program,  but  is  a  fundamental  limitation  of  estimation  capability 
in  the  presence  of  noise. 

No  derivatives  of  any  of  the  error  functions  to  be  extremized  are 
required  in  this  approach.  Instead,  direct  parabolic  interpolation  of  the 
appropriate  sampled  functions  is  employed  and  can  be  carried  through  several 
stages  to  the  desired  degree  of  accuracy  or  until  round-off  error  dominates. 
For  a  very  near  0  or  it,  the  approximation  of  8(a)  by  a  parabola  may  not  be 
adequate;  special  techniques  may  be  required  at  these  limits. 
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APPENDIX 

PROGRAM  FOR  ESTIMATION  OF  TONE  PARAMETERS 

Inputs  required  of  the  user  are  in 

line  10:  N,  number  of  data  points, 
line  20:  M,  size  of  FFT. 

The  program  is  configured  to  accept  up  to  N  =  8000  data  points  and  an  FFT  size 
up  to  M  =  16384.  The  user  can  also  change  the  weighting  from  flat  (in  line 
200)  to  whatever  weighting  is  of  interest.  The  appropriate  window  FFT  is 
undertaken  automatically,  by  means  of  lines  290  and  560.  The  initial  estimate 
of  a  and  the  plot  of  B(a)  are  completed  by  line  1010.  If  refined  estimates  of 
a  are  desired,  CONT  EXECUTE  must  be  performed,  and  can  be  repeated  for 
additional  refinement. 

The  terminology  DOUBLE  denotes  INTEGER  variables  in  BASIC  on  the  HP  9000 
computer.  Subroutine  SUB  B  computes  B(a)  and  coefficients  a  and  6  at  any 
frequency  a  of  interest.  Generation  of  data  for  the  examples  here  is 
accomplished  in  SUB  Data,  which  must  be  replaced  by  the  user  to  bring  in  his 
own  data. 
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10  H  =  2  5 

20  M =1024 

30  PRINT 

4Q  PRINT 

50  run  W<1 

60  RE DIM 

70  DOUBLE 

80  IF  M>2*N 

98  BEEP 

100  PRINT  "M 

110  PRUSE 

120  T  =  2 .  *  P  I  / M 

130  FOR  M s=0  TO  M/4 

140  Cos  <  Ms  )=COS  <  T*Ms  > 

150  NEXT  Ms 

160  MRT  X=(0.) 

170  MRT  Y  =  <  0  .  ) 

180  8=0. 

190  FOR  Ks= 1  TO  N 

200  W  k  =  1  . 

210  W  <  K  s  >  =  W  k 

220  S=S+Wk 

230  NEXT  Ks 

240  8=1. /S 

250  W 1 =W2=0 . 

260  FOR  Ks  = 1  TO  N 

270  T  =  W  <  K  s  >  *  S 

280  W  <  K  s  )  =  T 

290  Y  <  K  s  +  K  s  >  =  T 

300  T=T*Ks 

310  W 1  =  W 1  +  T 

320  W2=W2+T*Ks 

330  NEXT  Ks 

340  CRLL  D at a< N , Xd ( * ) > 

350  So=To=Se=Te=En=0 . 

360  FOR  K s  = 1  TO  N  STEP  2 

370  T 1 =  X d  < K  s  > 

380  T2=W  ( Ks ) *T 1 

390  X  <  Ks  >  =  T  2 

4  00  So  =  So  +  T  2 

410  To=To+T2*Ks 

420  En*En+Tl*T2 

430  NEXT  Ks 

440  FOR  K s  =  2  TO  N  STEP  2 

450  Tl=Xd<Ks> 

460  T2  =  W  <  Ks  >  *T 1 

470  X ( K  s ) =  T  2 

480  Se=Se+T2 

490  Te=Te+T2*Ks 

580  Er.=En  +  Tl*T2 

518  NEXT  Ks 


!  INTEGERS 

DECREASE  N." 

'  QUARTER-COSINE  TABLE 

i  SPECIFY  WEIGHTS,  k=l:N 

!  SCALE  WEIGHTS 
!  SO  THAT  SUM  =  1 

!  MOMENTS  OF 
!  WEIGHTS 

!  FILL  UP  DATA  ARRAY  Xd < 1 : N ) 


!  TOTAL  WEIGHTED  ENERGY 


!  NUMBER  OF  DATA  POINTS 
!  SIZE  OF  FFT ;  M  >  2N  REQUIRED 
"NUMBER  OF  DATA  POINTS  N  = " ; N 
"SIZE  OF  FFT  M  = " ; M 

000  >  ,  Xd  '•  1  :  8000  >  ,  X  <  0  :  1  6383  >  ,  Y  '•  0 :  16383 )  ,  Cos  <  0  :  4096  > 

W  <  1  :  N  >  ,  Xd  1  :  N ')  ,  X  <  0  :  M- 1  >  ,  Y  •;  0  :  M-  1  >  ,  Cos  <  0 :  M/4  > 

N, M, Ms, Ks, Js, M2 
THEN  120 


<=  2N 5  INCREASE  M  OR 
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520 

L  0  =  3  e  +  S  o 

!  MOMENTS 

5  3  O 

Ll=Te+To 

!  OF 

54  0 

L  0 1 =Se -So 

'  WEIGHTED 

550 

L  It  =  T  e - T o 

!  DATA 

560 

CALL  Fft  14>: 

M , Cos<*  > , X<*> , y< 

*  >  > 

570 

M  2  =  M  2 

5  8  0 

T  =  W  2  “  W  1  * W 1 

!  EVALUATE  B<a>  IN  X<0;  Mx; 

590 

X<0>=<L1*L1 

-2. *W1*L1*L0+W2* 

L  0  *  L  0  )  x  T 

6  0  0 

X  M2  =  <  L  1 1  *  L 1 1  -  2  .  *Wl*Llt*L  0 1  +  W2  *  L  0 1  *  L  0 1  >  /  T 

610 

FOR  M s  =  1  TO 

M2”  1 

620 

T 1 =X ( Ms > 

630 

T2  =  X  < M-Ms > 

640 

S  x  =  T  1  +  T  2 

650 

D  x  =  T  1  -  T  2 

660 

T  1  =  Y  <  M  s  > 

670 

T2  =  Y  < M-Ms  > 

630 

3  y  =  T  1  +  T  2 

6  90 

D  y  =  T  1  -  T  2 

700 

T  1  =  D  x  *  S  x  *  D  y 

710 

T  2  =  4 , - <  S  y  * S 

y+Dx*Dx > 

720 

X <  Ms >  =  < Sx*S 

x*(2. - S y  >  +  D  y *  D  y * 

2 .  +  S  y  )  +  T  1  +  T 1  )  /  T  2 

730 

NEXT  Ms 

740 

B i g  =  X C 0  ) 

!  SEARCH  FOR  MAXIMUM 

750 

JS  =  0 

760 

FOR  M  s  =  1  TO 

M2 

770 

T  =  X  <  M  s  ) 

730 

IF  T <  =  B i g  THEN  810 

790 

B  i  g  =  T 

!  MAXIMUM  VALUE  AND 

8  0  0 

Js=Ms 

1  LOCATION  IN  ARRAY 

8  1  0 

NEXT  Ms 

820 

IF  J s > 0  AND 

J s <  M 2  THEN  850 

330 

T  =  0 . 

340 

350 

GOTO  890 
T1=X< Js+1 > 

8  6  0 

T2=X< Js-1 > 

3  7  0 

T  =  . 5*<T1-T2 

>x<Big+Big-Tl-T2 

)  !  PARABOLIC  INTERPOLATION 

8  8  0 

B  i  g  =  B  i  g  +  .  2  5  *  <  T  1  -  T  2  >  *  T 

!  FOR  MAXIMUM  VALUE 

S  9  0 

ft  s  —  (.  J  s  +  T )  *  2 

.  *  P  I  •••'  M 

!  AND  LOCATION  OF  MAXIMUM 

9  0  0 

CALL  B < N , A s 

, W<*>, Xd<*> , fllph 

a ,  Be t-  a ,  B  =t> 

910 

PRINT  " INITIAL:  " ; " Js  =" ; Js 

; "  Big  =" ; Bi g; "  a  =" ; As; "  B C 

920 

G I  N I  T 

9  3  0 

PLOTTER  IS 

"GRAPHICS" 

940 

GRAPHICS  ON 

950 

WINDOW  0 . ,  M 

2, 0. 5  Ba 

96  0 

GRID  M 2 x 8 .  , 

B  ax  1 0 . 

970 

FOR  M s  =  0  TO 

M2 

9  8  0 

PLOT  MSpXCM 

s ) 

!  PLOT  B  (.  a  )  OVER  CO,  PI] 

9  9  0 

NEXT  Ms 

1000 

PENUP 

1010 

PAUSE 

1020 

GRAPHICS  OFF 
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1  0  3  0 
1  0  4  0 
1050 
1060 
1  070 
1030 
1  0  9  0 
1  1  00 
1110 
1  120 
1  130 
1  140 
1150 
1  160 
1  170 
1  180 
1  190 
1  200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1290 
1300 
1310 
1320 
1330 
1340 
1350 
1  3  6  0 
1  370 
1380 
1390 
1400 
1410 
1420 
1430 
1440 
1450 


De-l  t  a- 1  . /M  !  INITIAL  SEARCH  INCREMENT 

De  1  t  a=De  1  t  a*  .  1  !  FINE-GRAIN  SEARCH 

CALL  B  i  N  ,  A  s  -  D  e  1  t  a ,  N  *  >  ,  X  d  <  *  )  ,  A 1  p  h  a ,  B  e  t  a ,  B  am 
CALL  B  <  N  ,  As  +  De  lta,  U  <  *  >  ,  Xd  < *  >  j  A 1  pha,  Bet  a.  Bap  > 

T  = .  5  *  k  B  ap  -  B  am )  / B  a+  B  a  —  B  am  -  E  ap  ) 

IF  ABS(T)<1.  THEN  1100 

PRINT  "REFINED  INTERPOLATION  IS  BEYOND  EDGES  OF  SEARCH  INCREMENT:  ";T 

As=As+T*Del  t  a 

B  a=  B  a+  .  2  5  *  <  B  ap  -  B  am  >  *  T 

PRINT  "REFINED  a  =  " 5  As 

PRINT  "REFINED  E  <  a )  = " 5  B  a 

CALL  B < N , As, N < * ) , Xd < * ) , A 1 pha, Bet  a, Ea) 

PRINT  "MAXIMUM  B<a>  =";Ba 

Em i n=En-Ba  !  M I N I  MUM  ENERGY 

PRINT  A  1  pi  1  a  -  ;  A 1  pha  5  "  Bet  a  —  "  ;  Bet  a;  11  Em  i  n  =  "  ;  Em  i  n 

PRINT 
PAUSE 
GOTO  1040 
END 
! 

SUB  B ( DOUBLE  N , REAL  As , U < * > , Xd < * > , A1 pha, Beta, Ba) 

DOUBLE  Ks 

A2=As+As 

Ur  =  U i =  L  r  =  L i =0 . 

FOR  Ks=  1  TO  N 
T  u  =  U  < Ks ) 

Tx=Tw*Xd< K s  > 

T  1  =  A  s  *  K  s 
T  2= AL*Ks 

W  r  =  Ur  +  Tw*CQ$<!T2> 

Wi  =  Ui  -Tui*SINCT2> 

Lr=Lr+Tx*COS<T 1 > 

Li =Li -  T  x  *  S I N  <  T 1  > 

NEXT  Ks 
T 1  =  < 1 . -Ur  >  *Lr 
T  2  =  < 1 . +  U  r  )  *  L  i 
T  =  2 . / <  1  ■  - < H r  *  U r  +  N i *  W i > > 

A 1  p  h  a=  C  T 1  -  U  i  *  L  i  }  *  T 
Beta=c;Wi  *Lr-T2>*T 
B  a=  U  i  *  L  r  *  L  i 

B  a=  <.  T  1  *  L  r  +  T  2  *L  i  -  B  a-  B  a  >  *  T 
SUBEND 
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1460 

SUE  Ff t 14 (DOUBLE  N , F 

1470 

DOUBLE  N1 , N2, N3, N4, l 

1480 

DOUBLE  11,12,13,14, 

1490 

IF  N= 1  THEN  SUBEXIT 

1  500 

IF  N > 2  THEN  1580 

15  10 

R  =  X  0  >  +  X  <  1  > 

1520 

X  (  1  >  =  X  <  0  >  -  X  <  1  > 

1538 

X  (  0  )  =  H 

1540 

R  =  Y  <  O  >  +  Y  <  1  > 

1550 

Y< 1 >=Y<0)-Y< 1 > 

1560 

Y  <  0  >  =  R 

1570 

SUBEXIT 

1580 

N  1  =  N  /  4 

1590 

N  2  =  N 1  + 1 

1600 

N3=N2+1 

1610 

H  4  =  N  3  +  N  1 

1620 

L  o  g  2  n  =  1 . 4  4  2  7  *  L  O  G  <  N  > 

1  630 

FOR  11  =  1  TO  L o g 2 n 

1640 

I  2  =  2  A  <  L  o  g  2  n  -  I  1  > 

1650 

13=12+12 

1660 

I  4  =  N  /  I  3 

1670 

FOR  15=1  TO  12 

1680 

I  6  =  < 15-1 )* I 4+1 

1690 

IF  I  6 <  =  N 2  THEN  1730 

1700 

fll=-Cos(N4-I6-l> 

1710 

R  2  =  - C  o  s  < I6-N1-1 > 

1720 

GOTO  1750 

1730 

R  1  =  C  o s ( 16-1 > 

1746 

R2=“Cos<N3“I 6 - 1 > 

1750 

FOR  17=0  TO  N -13  ST 

1760 

18=17+15-1 

1770 

19=18+12 

1780 

T  1  =  X  <  1 8  > 

1790 

T  2  =  X  <  I  9  > 

1800 

T  3  =  Y  <  I  8  > 

1810 

T 4  =  Y <  ISO 

1820 

R  3  =  T 1-T2 

1830 

R4  =  T  3-T4 

1840 

X<  I8>  =  T1+T2 

1850 

Y  <  1 8  >  =  T3  +  T4 

I860 

X <  I9)=R1*R3-R2  +  R4 

1870 

Y<  I9)=R1+R4  +  R2  +  R3 

1  88  0 

NEXT  17 

1890 

NEXT  15 

1900 

NEXT  11 

1 6 ,  r 


*  >  ,  X  (  *  >  ,  Y  <  *  >  )  !  N  <  =  2  14=163  8  4  ;  0  S  U  E 
!  INTEGERS  <  2A31  =  2,147,483,64 
18, 19, 110, Ill, 112, 113, I  1 4 , L  <  0 :  13) 
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1910 

1 1  =  L  o  g  2  n  +  1 

1  920 

FOR  12=1  TO  14 

1930 

L  <  I  2  -  1  >  =  1 

1  940 

IF  I  2  >  L  o g 2  n  THEN  I960 

1  950 

L(  12-1  >=2M  11-12) 

1  9  6  0 

NEXT  12  * 

1970 

K  =  0 

1980 

FOR  11=1  TO  L(13) 

1  9  9  0 

FOR  12=11  TO  L C 1 2 )  STEP  L<13) 

2  0  0  0 

FOR  13=12  TO  L < 1 1 )  STEP  L(12) 

20  1  0 

FOR  14=13  TO  L  < 1 0  >  STEP  L < 1 1 ) 

2020 

FOR  15=14  TO  L <9 )  STEP  L< 10) 

2  0  3  O 

FOR  16=15  TO  L < 8 )  STEP  L < 9 > 

204  0 

FOR  17=16  T  i  J  L  (  7  )  STEP  L  <*  8  ) 

2050 

FOR  18=17  TO  L  < 6 )  STEP  L < 7 ) 

206  0 

FOR  19=13  TO  L  <  5 )  STEP  L  <  6  > 

2070 

FOR  110=19  TO  L  < 4  >  STEP  L < 5 ) 

2080 

FOR  111=110  TO  L<3>  STEP  L<4> 

2  0  9  0 

FOR  112=111  TO  L < 2 )  STEP  L<3) 

2100 

FOR  113=112  TO  L  < 1 )  STEP  L<2> 

2110 

FOR  114=113  TO  L  <  0 )  STEP  L<1> 

2120 

J  =  I  1  4  - 1 

2130 

IF  K  >  J  THEN  2200 

2140 

R  =  X  <  K  > 

2150 

X  <  K  >  =  X  <  J  ) 

2160 

X< J>=R 

2170 

R  =  V  <  K  > 

2130 

V  <  K  )  =  V  (  J  > 

2190 

Y  <  J  >  =  R 

2200 

K  =  K  + 1 

2210 

NEXT  114 

2228 

NEXT  113 

2230 

NEXT  112 

2240 

NEXT  Ill 

2250 

NEXT  IlO 

2260 

NEXT  19 

2270 

NEXT  18 

2230 

NEXT  17 

2290 

NEXT  16 

2300 

NEXT  15 

2310 

NEXT  14 

2320 

NEXT  13 

2330 

NEXT  12 

2340 

NEXT  I 1 

2350 

SUBEND 

2360 

i 

2370 

SUB  Dat a< DOUBLE  N, REfiL  Xd<*>> 

2380 

DOUBLE  Ks 

2390 

FOR  Ks= 1  TO  N 

2430 

Xd < Ks  > =COS  <  Ks )  + . 5*S I N  < Ks > 

2410 

NEXT  Ks 

2420 

SUBEND 

36 


TR  7785 


REFERENCES 

1.  L.  C.  Palmer,  "Coarse  Frequency  Estimation  Using  the  Discrete  Fourier 
Transform,"  IEEE  Trans,  on  Information  Theory,  vol.  IT-20,  Jan.  1974, 
pp.  104-109. 

2.  D.  C.  Rife,  "Digital  Tone  Parameter  Estimation  in  the  Presence  of 
Gaussian  Noise,"  Ph.D.  dissertation.  Polytechnic  Institute  of  Brooklyn, 
June  1973. 

3.  D.  C.  Rife  and  R.  R.  Boorstyn,  "Single  Tone  Parameter  Estimation  From 
Discrete-Time  Observations,"  IEEE  Trans,  on  Information  Theory, 

vol.  IT-20,  Sept.  1974,  pp.  591-598. 

4.  D.  C.  Rife  and  R.  R.  Boorstyn,  "Multiple  Tone  Parameter  Estimation  From 
Discrete  Time  Observations,"  BSTJ,  vol.  55,  no.  9,  Nov.  1976,  pp. 
1389-1410. 

5.  T.  Abatzoglou,  "A  Fast  Maximum  Likelihood  Algorithm  for  Frequency 
Estimation  of  a  Sinusoid  Based  on  Newton's  Method,"  IEEE  Trans,  on 
Acoustics,  Speech  and  Signal  Processing,  vol.  ASSP-33,  Feb.  1985,  pp. 
77-89. 

6.  E.  Plotkin  et  al.,  "Estimation  of  an  Initial  Phase  Over  Extremely  Short 
Record-Length  of  a  Sinewave  Signal,"  IEEE  Trans,  on  Instrumentation  and 
Measurement,  vol.  IM-34,  Dec.  1985,  pp.  624-629. 

7.  R.  Barrett  and  D.  McMahon,  "Comparison  of  Frequency  Estimators  for 
Underwater  Acoustic  Data,"  Journal  of  the  Acoustical  Society  of  America, 
vol.  79,  no.  5,  May  1986,  pp.  1461-1471. 

8.  R.  Kenefic,  "On  the  Estimation  of  an  Initial  Phase  Over  Extremely  Short 
Record-Length  of  a  Sinewave  Signal,"  submitted  to  the  IEEE  Trans,  on 
Instrumentation  and  Measurement,  August  1986. 


37 


TR  7785 


REFERENCES  (Cont'd) 

9.  R.  Kenefic  and  A.  H.  Nuttal 1 ,  "Maximum  Likelihood  Estimation  of  the 
Parameters  of  a  Tone  Using  Real  Discrete  Data,"  Journal  of  Oceanic 
Engineering,  vol.  OE-12,  no.  1,  Special  Issue  on  Underwater  Acoustic 
Signal  Processing,  January  1987. 

10.  A.  H.  Nuttall,  Accurate  Efficient  Evaluation  of  Cumulative  or  Exceedance 
Probability  Distributions  Directly  from  Characteristic  Functions,  NUSC 

Technical  Report  7023,  1  October  1983. 


38 


TR  7785 


INITIAL  DISTRIBUTION  LIST 

Addressee  No.  of  Copies 

ASN  ( RE&S)  1 
OUSDR&E  (Research  and  Advanced  Technology)  2 
DEPUTY  USDR8.E  (Res  &  Adv  Tech)  1 
DEPUTY  USDR&E  (Dir  Elect  &  Phys  Sc)  1 
ONR,  0NR-100,  -102,  -200,  -400,  -410,  411  (N.  Gerr) , 

-422,  -425AC ,  -430  9 
COMSPAWARSYSCOM,  SPAWAR  05  (W.  R.  Hunt)  2 
CNO ,  OP-098,  OP-941,  OP-951  3 
DIA  ( DT-2C)  10 
NRL ,  Code  5132  (Dr.  P.  B.  Abraham),  Code  5370 

(W.  Gabriel),  Code  5135  (N.  Yen),  (A.  A.  Gerlach)  5 
USRD  1 
NORDA,  Code  345  (R.  Wagstaff)  2 
USOC,  Code  240,  Code  241  2 
NAVSUBSUPACNLON  1 
NAVOCEANO,  Code  02  2 
NAVELECSYSCOM,  ELEX  03,  310  2 
NAVSEASYSCOM,  SEA-00,  -05R,  -06F,  63R,  -92R  5 
NAVAIRDEVCEN,  Warminster  1 
NAVAIRDEVCEN ,  Key  West  1 
NOSC ,  Code  8302,  Code  6565  (Library),  Code  713  (F.  Harris)  3 
NAVWPNSCEN  1 
NCSC ,  Code  724  1 
NAVCIVENGRLAB  1 
NAVSWC  1 
NAVSURFWPNCEN ,  Code  U31  1 
NISC,  1 
CNET,  Code  017  1 
CNTT  1 
NAVSUBSCOL  1 
NAVTRAEQUIPCENT,  Technical  Library  1 
NAVPGSCOL  1 
NAVWARCOL  1 
NETC  1 
APL/UW,  SEATTLE  1 
ARL/PENN  STATE,  STATE  COLLEGE  1 
CENTER  FOR  NAVAL  ANALYSES  (ACQUISITION  UNIT)  1 
DTIC  12 
DARPA,  A.  Ellinthorpe  2 
NOAA/ERL  1 
NATIONAL  RESEARCH  COUNCIL  1 
WOODS  HOLE  OCEANOGRAPHIC  INSTITUTION  (Dr.  R.  C.  Spindel)  2 
ENGINEERING  SOCIETIES  LIB,  UNITED  ENGRG  CTR  1 
NATIONAL  INSTITUTE  OF  HEALTH  1 
ARL,  UNI V  OF  TEXAS  1 


TR  7785 


INITIAL  DISTRIBUTION  LIST  (Cont'd) 

Addressee  No.  of  Copies 

MARINE  PHYSICAL  LAB,  SCRIPPS  1 
UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO  1 
NAVSURWEACTR  1 
DELSI  1 
DIRECTOR  SACLANT  ASW  RES  CEN  1 
COM  SPACE  &  NAV  WAR  SYS  COM  1 
COM  NAVAL  PERSONNEL  R&D  CENTER  1 
COM  NAV  SUB  COLLEGE  1 
B-K  DYN  INC  1 
BBN,  Arlington,  V A  (Dr.  H.  Cox)  1 
BBN,  Cambridge,  MA  (H.  Gish)  1 
BBN,  New  London,  CT  (Dr.  P.  Cable)  1 
EWASCTRI  1 
MAR,  INC,  East  Lyme,  CT  1 
HYDROINC  (D.  Clark)  1 
SUMRESCR  (M.  Henry)  1 
ANALTECHNS,  N.  Stonington,  CT  (H.  Jarvis)  1 
ANALTECHNS,  New  London,  CT  1 
EDOCORP  (J.  Vincenzo)  1 
TRA  CORP.,  Austin,  TX  (Dr.  T.  Leih,  J.  Wilkinson)  2 
TRA  CORP.,  Groton,  CT  1 
NETS  (R.  Medeiros)  1 
GESY,  D.  Bates  1 
SONALYSTS,  Waterford,  CT  (J.  Morris)  1 
ORI  CO,  INC.  (G.  Assard)  1 
HUGHES  AIRCRAFT  CO.  (S.  Autrey)  1 
MIT  (Prof.  A.  Baggaroer)  1 
RAYTHEON  CO.  (J.  Bartram)  1 
Dr.  Julius  Bendat,  833  Moraga  Dr,  Los  Angeles,  CA  1 
COOLEY  LABORATORY  (Prof.  T.  Birdsall)  1 
PROMETHEUS,  INC.  (Dr.  J.  S.  Byrnes)  1 
ROYAL  MILITARY  COLLEGE  OF  CANADA  (Prof.  Y.  T.  Chan)  1 
UNIV.  OF  FLORIDA  (D.  C.  Childers)  1 
SANDIA  NATIONAL  LABORATORY  (J.  Claasen)  1 
COGENT  SYSTEMS,  INC.  (J.  P.  Costas)  1 
IBM  CORP.  (G.  Demuth)  1 
UNIV.  OF  STRATHCLYDE,  CLASGOW,  SCOTLAND  (Prof.  T.  Durrani)  1 
ROCKWELL  INTERNATIONAL  CORP.  (L.  T.  Einstein  and 

Dr.  D.  F.  Elliott)  2 
GENERAL  ELECTRIC  CO.  (Dr.  M.  Fitelson)  1 
HONEYWELL,  INC.  (D.  M.  Goodfellow)  1 
UNIV.  OF  TECHNOLOGY,  LOUGHBOROUGH,  LEICESTERSHIRE, 

ENGLAND  (Prof.  J.  W.  R.  Griffiths)  1 
HARRIS  SCIENTIFIC  SERVICES  (B.  Harris)  1 
UNIV  OF  CALIFORNIA,  SAN  DIEGO  (Prof.  C.  W.  Helstrom)  1 
EG&G  (Dr.  J.  Hughen)  1 


TR  7785 


INITIAL  DISTRIBUTION  LIST  (Cont'd) 

Addressee  No.  of  Copies 

BELL  COMMUNICATIONS  RESEARCH  (J.  F.  Kaiser)  1 
UNIV.  OF  RI  (Prof.  S.  Kay,  Prof.  L.  Scharf, 

Prof.  D.  Tufts)  3 
DREXEL  UNIV.  (Prof.  S.  Kesler)  1 
UNIV.  OF  CT  (Prof.  C.  H.  Knapp)  1 
APPLIED  SEISMIC  GROUP  (R.  Lacoss)  1 
ADMIRALTY  RESEARCH  ESTABLISHMENT,  ENGLAND 

(Dr.  L.  J.  Lloyd)  1 
NAVAL  SYSTEMS  DIV.,  SIMRAD  SUBSEA  A/S,  NORWAY  (E.  B.  Lunde)  1 
DEFENCE  RESEARCH  ESTAB.  ATLANTIC,  DARTMOUTH,  NOVA  SCOTIA 

CANADA  (B.  E.  Mackey,  Library)  1 
MARTIN  MARIETTA  AEROSPACE  (S.  L.  Marple)  1 
PSI  MARINE  SCIENCES  (Dr.  R.  Mellen)  1 
Dr.  D.  Middleton,  127  E.  91st  St.  NY,  NY  10128  1 
Dr.  P.  Mikhalevsky,  803  W.  Broad  St. 

Falls  Church,  VA  22046  1 
CANBERRA  COLLEGE  OF  ADV.  EDUC.,  AUSTRALIA  2616  (P.  Morgan)  1 
NORTHEASTERN  UNIV.,  (Prof.  C.  L.  Nikias)  1 
ASTRON  RESEARCH  &  ENGR,  (Dr.  A.  G.  Piersol)  1 
WESTINGHOUSE  ELEC.  CORP.  (Dr.  H.  L.  Price)  1 
M/A-COM  GOVT  SYSTEMS,  (Dr.  R.  Price)  1 
DALHOUSIE  UNIV.  HALIFAX,  NOVA  SCOTIA  (Dr.  B.  Ruddick)  1 
NATO  SACLANT  ASW  RESEARCH  CENTER,  LIBRARY,  APO  NY,  NY  09019  1 
YALE  UNIV.  (Prof.  P.  M.  Schultheiss)  1 
NATIONAL  RADIO  ASTRONOMY  OBSERVATORY  (F.  Schwab)  1 
DEFENSE  SYSTEMS,  INC  (Dr.  G.  S.  Sebestyen)  1 
NATO  SACLANT  ASW  RESEARCH  CENTRE  (Dr.  E.  J.  Sullivan)  1 
PENN  STATE  UNIV,  APPLIED  RESEARCH  LAB.  (F.  W.  Symons)  1 
NAVAL  PG  SCHOOL,  (Prof.  C.  W.  Therrien)  1 
DEFENCE  RESEARCH  ESTABLISHMENT  PACIFIC,  VICTORIA,  B.C. 

CANADA  VOS  1  BO  (Dr.  D.  J .  Thomson)  1 
Robert  J.  Urick,  11701  Berwick  Rd,  Silver  Spring,  MD  20904  1 
RCA  CORP  (H.  Urkowitz)  1 
USEA  S.P.A.  LA  SPEZIA,  ITALY  (H.  Van  Asselt)  1 
TEL-AVIV  UNIV.  ISRAEL  (Prof.  E.  Weinstein)  1 
COAST  GUARD  ACADEMY  (Prof.  3.  3.  Wolcin)  1 
SPACE  PHYSICS  LAB,  UNIV  OF  ALBERTA,  EDMONTON,  CANADA 

(K.  L.  Yeung)  1 
COLUMBIA  RESEARCH  CORP.  (W.  Hahn)  1 
SAIC ,  NLON  (Dr.  F.  DiNapoli)  1 


U2267 12 


